//===================================================== file = genemp.c =====
//=  Program to generate empirically distributed random variables           =
//===========================================================================
//=  Notes: 1) Writes to a user specified output file                       =
//=         2) Generates user specified number of samples                   =
//=         3) The empirical distribution is assume to exist in a file      =
//=            named "dist.dat".  The format of dist.dat is set of doubles  =
//=            in column order.  Each double is the probability of          =
//=            occurence and the value of the RV.  The probability          =
//=            values must sum to 1.0.  See the example below.              =
//=-------------------------------------------------------------------------=
//= Example user input:                                                     =
//=                                                                         =
//=   ----------------------------------------- genemp.c -----              =
//=   -  Program to generate empirically distributed random  -              =
//=   -  variables (empirical distribution in file dist.dat) -              =
//=   --------------------------------------------------------              =
//=   Output file name ===================================> output.dat      =
//=   Random number seed =================================> 1               =
//=   Number of samples to generate ======================> 5               =
//=   --------------------------------------------------------              =
//=   -  Building the empirical CDF                          -              =
//=   --------------------------------------------------------              =
//=   --------------------------------------------------------              =
//=   -  Generating samples to file                          -              =
//=   --------------------------------------------------------              =
//=   --------------------------------------------------------              =
//=   -  Done!                                                              =
//=   --------------------------------------------------------              =
//=-------------------------------------------------------------------------=
//= Example input file (this is the required file "dist.dat"                =
//=                                                                         =
//=   0.25  5.00                                                            =
//=   0.45  2.00                                                            =
//=   0.15  1.00                                                            =
//=   0.10  0.50                                                            =
//=   0.05  0.25                                                            =
//=                                                                         =
//= Example output (from above input file and user input):                  =
//=                                                                         =
//=   5.000000                                                              =
//=   5.000000                                                              =
//=   1.000000                                                              =
//=   2.000000                                                              =
//=   2.000000                                                              =
//=-------------------------------------------------------------------------=
//=  Build: gcc genemp.c, bcc32 genemp.c, cl genemp.c                       =
//=-------------------------------------------------------------------------=
//=  Execute: genexp (must have a file "dist.dat" in same directory)        =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (12/30/00) - Genesis (from genexp.c)                      =
//=           KJC (01/29/01) - Changed Num_entries to be an int variable    =
//=           KJC (03/12/03) - Added Jain's RNG for finer granularity       =
//===========================================================================

//----- Include files -------------------------------------------------------
#include <stdio.h>              // Needed for printf()
#include <stdlib.h>             // Needed for exit() and ato*()

//----- Defines -------------------------------------------------------------
#define DIST_FILE  "dist.dat"   // Distribution file
#define MAX_ENTRY       1000    // Maximum number of lines in DIST_FILE

//----- Globals -------------------------------------------------------------
struct rv                       // Structure for an empirical RV
{
    double cdf_val;               // *** Probability value
    double rv_val;                // *** RV value
};
struct rv CDF[MAX_ENTRY];       // The CDF built from DIST_FILE
int       Num_entries;          // The number of entries in the CDF

//----- Function prototypes -------------------------------------------------
double emp(void);               // Returns an empirical random variable
double rand_val(int seed);      // Jain's RNG

//===== Main program ========================================================
void main(void)
{
    FILE     *fp_dist;              // File pointer to distribution file
    FILE     *fp_out;               // File pointer to output file
    char     instring1[80];         // Input string #1
    char     instring2[80];         // Input string #2
    double   num_samples;           // Number of samples to generate
    double   emp_rv;                // Exponential random variable
    int      i;                     // Loop counter

    // Output banner
    printf("----------------------------------------- genemp.c ----- \n");
    printf("-  Program to generate empirically distributed random  - \n");
    printf("-  variables (empirical distribution in file dist.dat) - \n");
    printf("-------------------------------------------------------- \n");

    // Open the distribution file
    fp_dist = fopen(DIST_FILE, "r");
    if (fp_dist == NULL)
    {
        printf("ERROR in opening the distributin file (%s) \n", DIST_FILE);
        exit(1);
    }

    // Build the CDF from the distribution file
    Num_entries = 0;
    fscanf(fp_dist, "%s %s \n", instring1, instring2);
    CDF[Num_entries].cdf_val = atof(instring1);
    CDF[Num_entries].rv_val = atof(instring2);
    while (!feof(fp_dist))
    {
        Num_entries++;
        fscanf(fp_dist, "%s %s \n", instring1, instring2);
        CDF[Num_entries].cdf_val = atof(instring1) + CDF[Num_entries - 1].cdf_val;
        CDF[Num_entries].rv_val = atof(instring2);
    }

    // Verify that last cdf_val in the CDF is 1.0 (sum of probabilities = 1.0)
    if (CDF[Num_entries].cdf_val != 1.0)
    {
        printf("ERROR - sum of probabilities in %s is %f (must be 1.0) \n",
               DIST_FILE, CDF[Num_entries].cdf_val );
        exit(1);
    }

    // Prompt for output filename and then create/open the file
    printf("Output file name ===================================> ");
    scanf("%s", instring1);
    fp_out = fopen(instring1, "w");
    if (fp_out == NULL)
    {
        printf("ERROR in creating output file (%s) \n", instring1);
        exit(1);
    }

    // Prompt for random number seed and then use it
    printf("Random number seed =================================> ");
    scanf("%s", instring1);
    rand_val((int) atoi(instring1));

    // Prompt for number of samples to generate
    printf("Number of samples to generate ======================> ");
    scanf("%s", instring1);
    num_samples = atoi(instring1);

    // Output message and build the empriical CDF
    printf("-------------------------------------------------------- \n");
    printf("-  Building the empirical CDF                          - \n");
    printf("-------------------------------------------------------- \n");

    // Output message and generate interarrival times
    printf("-------------------------------------------------------- \n");
    printf("-  Generating samples to file                          - \n");
    printf("-------------------------------------------------------- \n");
    for (i = 0; i < num_samples; i++)
    {
        emp_rv = emp();
        fprintf(fp_out, "%f \n", emp_rv);
    }

    // Output message and close the distribution and output files
    printf("-------------------------------------------------------- \n");
    printf("-  Done! \n");
    printf("-------------------------------------------------------- \n");
    fclose(fp_dist);
    fclose(fp_out);
}

//===========================================================================
//=  Function to generate empirically distributed random variables          =
//=    - Input: None                                                        =
//=    - Output: Returns with empirically distributed random variable       =
//===========================================================================
double emp(void)
{
    double z;                     // Uniform random number (0 < z < 1)
    double emp_value;             // Computed exponential value to be returned
    int    i;                     // Loop counter

    // Pull a uniform random number (0 < z < 1)
    z = rand_val(0);

    // Map z to empirical distribution and get empirical RV value
    for (i = 0; i <= Num_entries; i++)
        if (z <= CDF[i].cdf_val)
        {
            emp_value = CDF[i].rv_val;
            break;
        }

    // Return the empirical RV value
    return(emp_value);
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=========================================================================
double rand_val(int seed)
{
    const long  a =      16807;  // Multiplier
    const long  m = 2147483647;  // Modulus
    const long  q =     127773;  // m div a
    const long  r =       2836;  // m mod a
    static long x;               // Random int value
    long        x_div_q;         // x divided by q
    long        x_mod_q;         // x modulo q
    long        x_new;           // New x value

    // Set the seed if argument is non-zero and then return zero
    if (seed > 0)
    {
        x = seed;
        return(0.0);
    }

    // RNG using integer arithmetic
    x_div_q = x / q;
    x_mod_q = x % q;
    x_new = (a * x_mod_q) - (r * x_div_q);
    if (x_new > 0)
        x = x_new;
    else
        x = x_new + m;

    // Return a random value between 0.0 and 1.0
    return((double) x / m);
}
